A 3D Elastoplastic Constitutive Model Considering Progressive Damage Behavior for Thermoplastic Composites of T700/PEEK

Due to their excellent mechanical properties, the carbon fiber-reinforced polymer composites (CFRPs) of thermoplastic resins are widely used, and an accurate constitutive model plays a pivotal role in structural design and service safety. A two-parameter three-dimensional (3D) plastic potential was obtained by considering both the deviatoric deformation and the dilatation deformation associated with hydrostatic stress. The Langmuir function was first adopted to model the plastic hardening behavior of composites. The two-parameter 3D plastic potential, connected to the Langmuir function of plastic hardening, was thus proposed to model the constitutive behavior of the CFRPs of thermoplastic resins. Also, T700/PEEK specimens with different off-axis angles were subjected to tensile loading to obtain the corresponding fracture surface angles of specimens and the load–displacement curves. The two unknown plastic parameters in the proposed 3D plastic potential were obtained by using the quasi-Newton algorithm programmed in MATLAB, and the unknown hardening parameters in the Langmuir function were determined by fitting the effective stress-plastic strain curve in different off-axis angles. Meanwhile, the user material subroutine VUMAT, following the proposed constitutive model, was developed in terms of the maximum stress criterion for fiber failure and the LaRC05 criterion for matrix failure to simulate the 3D elastoplastic damage behavior of T700/PEEK. Finally, comparisons between the experimental tests and the numerical analysis were made, and a fairly good agreement was found, which validated the correctness of the proposed constitutive model in this work.


Introduction
Carbon fiber-reinforced polymer composites (CFRPs) have been integrated into many industrial applications for their structural parts with the requirement of a lightweight design and high reliability.Thermosetting resins, such as epoxy and others, have been widely used in various industry sectors.The thermoplastic resins, especially the polyetheretherketone (PEEK), have a higher impact resistance, better fracture toughness, and the advantage of waste recycling compared to the thermosetting resin [1][2][3].Thus, the CFRPs of thermoplastic resins have great potential to be widely used in future aerospace, automotive, and civilian products.The stress-strain relationships of CFRPs with thermoplastic resins tend to be nonlinear and sensitive to the plastic strain [4,5], herein an accurate constitutive model is one of the most critical considerations in its structural integrity.
Numerous studies on constitutive behavior have been conducted on the thermosetting composites, whereas relatively few studies have been conducted on the thermoplastic composites.The research methods of the former can provide available guidance for the latter.It is widely accepted that the plastic potential plays a key role in the constitutive model due to its calculation of the effective stress of materials [6,7].Based on a simple flow rule and the Hill plastic potential, Sun and Chen [8] first proposed a plastic potential by considering a state of plane stress.The plastic potential includes one undetermined parameter to characterize the level of plastic deformation developed under shear loading compared to transverse loading, and the single-parameter 2D plasticity model was employed to describe the nonlinear behavior of fibrous composites [9][10][11].Sun and Chen [12] assumed afterward that the composites were transversely isotropic, and a 3D plastic potential was established with two undetermined parameters.It is always desirable to have fewer parameters in the constitutive model by comparing the two-parameter 3D plastic potential with the strain energy function for transversely isotropic elastic solids.Thus, two undetermined parameters in the 3D plastic potential were reduced to one by Weeks and Sun [13], and a single-parameter 3D model was established.Cho et al. [14] introduced a 3D plastic potential, combining the functions related to dilatational deformations, and a simplified 2D plastic potential with two unknown parameters was given to perform the numerical prediction.Also, the validating experiment of unidirectional carbon/epoxy (AS4/3501-6) was carried out in this study.Chen and Suo [15] constructed an elastoplastic constitutive model to describe the nonlinear compressive stress-strain behavior of the fiber-reinforced polymer composites (FRPs) of T800H/3631 and HTS40/PA6, and the model was developed under the single-parameter 2D plastic potential.T700/PPEK is a high-performance CFRP of thermoplastic resins.Yao et al. [16] carried out a numerical and experimental study to investigate the effect of interference fit on the bearing strength of the T700/PEEK riveted joint.The single-parameter 2D plastic potential was chosen, and the unknown single parameter was obtained using the off-axis test.No works were seen in the open literature on combining the single-parameter 3D model with the functions related to the dilatational deformations for the simulation of the elastoplastic behavior for the CFRPs, and no similar works were seen for the CFRPs of thermoplastic resins.
The hardening laws were utilized to express the master effective stress-plastic strain relation, and the unknown coefficients in the hardening laws can be obtained by performing the off-axis tensile tests [17,18].Ding et al. [19] performed a 3D elastoplastic analysis of the interlaminar stresses for the AS4/PEEK composite laminate by using a single-parameter 3D plastic potential, and the master effective stress-plastic strain relation was fitted with a power-law function.An elastoplastic damage model considering the cohesive matrix interface layers was proposed for the composite laminates of AS4/PEEK by Mandal and Chakrabarti [20].It was pointed out that a hardening was considered in the shear direction, and there were two shear hardening parameters given in this study.Chen and Suo [15] concluded that the effective plastic strain could be expressed by the effective stress, initial elastic response stress, and two hardening parameters, which assume that the hardening behavior obeys a power law with an index independent of fiber orientation, tension, and compression.However, further improvements can be made by changing the functional form to enhance the accuracy of the numerical prediction.Vyas et al. [21] developed an exponential strain-hardening model with an explicit integration finite element code and studied the nonlinear mechanical response of AS4/55A under different transverse loads, showing a good agreement between the simulation and experimental results.A pressuredependent elastoplastic constitutive model was established for the unidirectional CFRP laminates of IM7/8552 and AS4/PEEK by Ren Rui et al. [22].The hardening law was expressed by a piecewise function composed of four terms, including three exponential functions, the curves of predictions agreed well with the tests, but there were seven unknown coefficients to be determined in this law.These functional forms of hardening law in the aforementioned literature were generally obtained using the simple linear combinations of exponential functions.The relation of the dependent variable-independent variable in the Langmuir isothermal adsorption function [23][24][25] shows similar patterns to that of the effective stress-effective plastic strain in the constitutive model.In addition, there is a lack of research on the application of the Langmuir function in the hardening law for the CFRPs of both thermosetting and thermoplastic resins.
Finite Element Analysis (FEA) plays a crucial role in the development of composite materials, which have been widely embraced by the industry and academic research [26,27].Chen et al. [6] developed a user-defined subroutine within ABAQUS/Standard to analyze the plastic behavior of composite materials.This subroutine incorporates a single-parameter plastic potential, Hashin failure criteria, and exponential damage evolution to predict the plastic behavior, damage initiation, and failure process of AS4/PEEK.Liu et al. [18] examined the three-point bending behavior of AS4/PEEK through a comprehensive approach involving both experimental and numerical methods.In addition to conducting experimental research, a numerical model was constructed to enhance the understanding of damage mechanisms and evolution.Din et al. [28] developed a UMAT subroutine in ABAQUS/Standard finite element software to investigate the mechanical properties of open-hole laminates, considering elastic-plastic continuous damage.The study utilized a one-parameter flow rule for orthotropic plasticity to capture the nonlinear behavior of the fiber-reinforced composite with a thermoplastic matrix.The Puck theory was employed as the failure criterion, along with an exponential damage evolution rule to predict damage onset and evolution.The analysis results demonstrate a notable enhancement in predicting the mechanical behavior of continuous unidirectional composite materials.Harpreet et al. [29] proposed an elastoplastic damage model for 3D fiber-reinforced plastic (FRP) composites to simulate progressive damage and induced inelastic deformation resulting from low-velocity impacts.The model utilized 3D plastic potential for plastic surface growth and implemented an exponential softening model for predicting damage growth.Additionally, delamination behavior was considered by incorporating surface-based cohesive interactions between adjacent plies.The results demonstrated a close match between the predicted behavior and the experimental observations.In recent years, numerous failure criteria have been proposed, with the LaRC05 theory being recognized as one of the most accurate based on the second World-Wide Failure Exercise (WWFE-II).Wang et al. [30] developed a progressive damage model using the 3D LaRC failure criterion in conjunction with cohesive elements.They utilized an energy-based damage evolution method to forecast the behavior of single-lap thin-ply laminated composite-bolted joints under tension loading.Ma et al. [31] conducted a series of compression tests using specimens at different off-axis angles to study the damage initiation and progression mechanisms in unidirectional AS4/PEEK composites.The off-axis compression failure envelope, evaluated based on the LaRC05 criteria, was compared to the experimental results, demonstrating the ability of the LaRC05 criteria to predict failure accurately.
To overcome the aforementioned issues and limitations, the CFRPs of thermoplastic resins-T700/PEEK-were regarded as the research object.A two-parameter 3D plastic potential was given by combining the single-parameter 3D plastic potential with the function terms related to the dilatational deformations, induced from the hydrostatic pressure.The Langmuir function was first chosen to express the hardening law of CFRPs.Furthermore, off-axis tensile tests were conducted to obtain the unknown parameters both in the plastic potential and the hardening law.Finally, the prediction of numerical simulations, considering the 3D elastoplastic damage behavior, was compared to the results of the off-axis test to validate the availability of the proposed constitutive model.

Material and Specimen
The unidirectional thermoplastic composites used in this study was T700/PEEK thermoplastic prepreg produced by Ningbo Material Technology and Engineering Research Institute, Ningbo, Zhejiang.It was fabricated by lying prepregs in a mental mold and then curing in a vacuum sulphuration at 380 • C and 3 MPa, as shown in Figure 1.There are seven samples with different fiber layups to be manufactured, namely [0] 8 , [15] 20 , [30] 16 .After cooling to room temperature, the specimens were cut to have off-axis angles, with respect to the fiber direction.The test specimens were cut by using a five-axis linkage ultra-high pressure water jet cutting machine at 27 MPa and 800 mm/min to ensure that the samples did not have serious cutting damage or interlayer delamination.The size of test specimens for the [0] 8 layup is 250 mm × 15 mm × 1 mm, the size of test specimens for the [90] 16 layup is 175 mm × 25 mm × 2 mm, and the size of test specimens for the other layups is 250 mm × 25 mm × 2.5 mm.Reinforcement plates are attached to both ends of the test specimen, the size of reinforcement plates for the [0] 8 layup is 56 mm × 15 mm × 2 mm, and the size of reinforcement plates for the other layups is 25 mm × 25 mm × 2 mm.At least five specimens were tested for each sample.
Materials 2024, 17, x FOR PEER REVIEW using a five-axis linkage ultra-high pressure water jet cutting machine at 27 MPa a mm/min to ensure that the samples did not have serious cutting damage or int delamination.The size of test specimens for the [0]8 layup is 250 mm × 15 mm × 1 m size of test specimens for the [90]16 layup is 175 mm × 25 mm × 2 mm, and the size specimens for the other layups is 250 mm × 25 mm × 2.5 mm.Reinforcement pla attached to both ends of the test specimen, the size of reinforcement plates for layup is 56 mm × 15 mm × 2 mm, and the size of reinforcement plates for the other is 25 mm × 25 mm × 2 mm.At least five specimens were tested for each sample.

Off-Axis Testing and Results
The experiment was conducted using the universal testing machine (MTS 37 Frame; MTS Systems Corp., Eden Prairie, MN, USA), as shown in Figure 2, and it tor force capacity was 100 kN.An extensometer with a 50 mm gauge length was eq with the machine for specimen tests.Off-axis tensile tests were carried out accor the standard ASTM D3039 [32] and the loading was applied under displacement at a rate of 2 mm/min.The frequency of data acquisition was set at 20 Hz.

Off-Axis Testing and Results
The experiment was conducted using the universal testing machine (MTS 370 Load Frame; MTS Systems Corp., Eden Prairie, MN, USA), as shown in Figure 2, and its actuator force capacity was 100 kN.An extensometer with a 50 mm gauge length was equipped with the machine for specimen tests.Off-axis tensile tests were carried out according to the standard ASTM D3039 [32] and the loading was applied under displacement control at a rate of 2 mm/min.The frequency of data acquisition was set at 20 Hz.The experimental data was processed according to ASTM D3039.The calculation formula for tensile stress is given by: where σ i denotes the tensile stress of the ith data point, Pi is the load of ith data point, and A represents the average cross-sectional area of the specimen.
The formula for calculating the elastic modulus is: in which E chord refers to the tensile chord modulus of elasticity; Δσ is the deviation value of applied tensile stress between two strain points on the extensometer; and Δε represents the deviation value of strain between two strain points.
Let the direction of uniaxial tensile loading be the x-axis, which is parallel to the longitudinal direction of rectangular specimens, and the average experimental axial stressstrain (σx − εx) curves of the off-axis tensile test are shown in Figure 3. From the stressstrain curves, we can see that when the off-axis angle is 0°, the mechanical behavior of the specimens is linearly elastic before its failure.As the off-axis angle gradually increases, the nonlinear behavior and ultimate fracture strain of specimens both show a decreasing trend.When the off-axis angle ranges from 60° to 90°, the mechanical behavior of the specimens only exhibits a small amount of approximately linearly elastic behavior, and then fracture failure occurs without apparent plastic deformation.This is because when the offaxis angle is large, the fiber cannot bear the majority of applied loads, resulting in the fracture failure of the matrix prior to the exhibition of plastic deformation.The experimental data was processed according to ASTM D3039.The calculation formula for tensile stress is given by: where σ i denotes the tensile stress of the ith data point, P i is the load of ith data point, and A represents the average cross-sectional area of the specimen.The formula for calculating the elastic modulus is: in which E chord refers to the tensile chord modulus of elasticity; ∆σ is the deviation value of applied tensile stress between two strain points on the extensometer; and ∆ε represents the deviation value of strain between two strain points.
Let the direction of uniaxial tensile loading be the x-axis, which is parallel to the longitudinal direction of rectangular specimens, and the average experimental axial stressstrain (σ x − ε x ) curves of the off-axis tensile test are shown in Figure 3. From the stressstrain curves, we can see that when the off-axis angle is 0 • , the mechanical behavior of the specimens is linearly elastic before its failure.As the off-axis angle gradually increases, the nonlinear behavior and ultimate fracture strain of specimens both show a decreasing trend.When the off-axis angle ranges from 60 • to 90 • , the mechanical behavior of the specimens only exhibits a small amount of approximately linearly elastic behavior, and then fracture failure occurs without apparent plastic deformation.This is because when the off-axis angle is large, the fiber cannot bear the majority of applied loads, resulting in the fracture failure of the matrix prior to the exhibition of plastic deformation.
The fracture surface angles are shown for the test specimens with different layups in Table 1.Due to the tensile failure mode of longitudinal splitting and fiber breakage in sample [0] 8 , the failure angle is set as zero.In terms of the off-axis angles among different specimens, the fracture surface angles increase accordingly, approximately the same value as the off-axis angle of the sample.The tested specimens with off-axis tension fracture failure are shown in Figure 4.
specimens is linearly elastic before its failure.As the off-axis angle gradually increases, the nonlinear behavior and ultimate fracture strain of specimens both show a decreasing trend.When the off-axis angle ranges from 60° to 90°, the mechanical behavior of the specimens only exhibits a small amount of approximately linearly elastic behavior, and then fracture failure occurs without apparent plastic deformation.This is because when the offaxis angle is large, the fiber cannot bear the majority of applied loads, resulting in the fracture failure of the matrix prior to the exhibition of plastic deformation.The fracture surface angles are shown for the test specimens with different layups in Table 1.Due to the tensile failure mode of longitudinal splitting and fiber breakage in sample [0]8, the failure angle is set as zero.In terms of the off-axis angles among different specimens, the fracture surface angles increase accordingly, approximately the same value as the off-axis angle of the sample.The tested specimens with off-axis tension fracture failure are shown in Figure 4.

Constitutive Model Coupling with 3D Elastoplastic Damage
From the results of the off-axis tensile tests, it can be seen that the CFRPs of T700/PEEK laminates exhibit a significant nonlinear mechanical response and typical tensile failure mode due to plastic deformation under the off-axis tensile load.In order to simulate the above mechanical behavior accurately, a 3D elastoplastic constitutive model of the CFRPs of thermoplastic resins was established by considering the elastic deformation, plastic deformation, damage initiation, and damage evolution.

Elasticity Description
By assumption of the transverse isotropy for unidirectional composite laminates, the expression of stress σ during the linear elastic stage is: where C0 represents the stiffness matrix of the initial undamaged composite laminates; ε e denotes the elastic strain tensor, and Equation ( 3) can be specifically given by:

Constitutive Model Coupling with 3D Elastoplastic Damage
From the results of the off-axis tensile tests, it can be seen that the CFRPs of T700/PEEK laminates exhibit a significant nonlinear mechanical response and typical tensile failure mode due to plastic deformation under the off-axis tensile load.In order to simulate the above mechanical behavior accurately, a 3D elastoplastic constitutive model of the CFRPs of thermoplastic resins was established by considering the elastic deformation, plastic deformation, damage initiation, and damage evolution.

Elasticity Description
By assumption of the transverse isotropy for unidirectional composite laminates, the expression of stress σ during the linear elastic stage is: where C 0 represents the stiffness matrix of the initial undamaged composite laminates; ε e denotes the elastic strain tensor, and Equation ( 3) can be specifically given by: and where the parameters E ij (i, j = 1, 2, 3, i = j) are the Young's modulus; ν ij (i, j = 1, 2, 3, i ̸ = j) are the Poisson 's ratios; G ij (i, j = 1, 2, 3, i ̸ = j) are the shear modulus.σ ij (i, j = 1, 2, 3) and ε e ij (i, j = 1, 2, 3) are the stress and strain components, respectively.

Plasticity Description
For the test specimens, the applied external force is resisted by the effective and undamaged area of the material.Therefore, it can be reasonably assumed that plastic deformation occurs in the undamaged region of the material.
The plastic yield criterion proposed by Sun and Chen [8,12] can effectively simulate the plastic mechanical behavior of the composite laminates under quasi-static conditions.The form of this model is adopted: where F p (σ) refers to the plastic potential; σ(ε p ) represents the hardening function.
Weeks et al. [13] assumed that the composite is transversely isotropic and only linearly elastic deformation occurred in its fiber direction.Also, the effect of dilatational deformation on the plastic deformation was considered by Cho et al. [14].By combining the single-parameter 3D plastic potential [13] and the dilatational function terms [14].A two-parameter 3D plastic potential is then proposed: where the values of a and b are both related to the properties of the composite material.As the dilatational deformations are ignored, namely b = 0, the proposed two-parameter 3D plastic potential in Equation ( 8) reduces to a single-parameter 3D plastic potential [13,33].Defining effective stress as: Assuming that the plastic deformation of the materials satisfies the associated flow rule, the incremental component of plastic strain can be expressed by: where dλ is the non-negative plastic multiplier of the entire plastic loading history; the gradient vector Combining Equations ( 8) and (10), we obtain: The hardening law is utilized to describe the yield stress evolution of the materials after entering the plastic phase.Sun and Chen [8,12] derived the effective plastic strain increment by defining the increment of effective plastic work per unit volume, dW p .They assumed that the product of effective stress and the increment of effective plastic strain is equal to the sum of the products of stress components in each direction and its corresponding increment of plastic strain.Finally, the relations can be given by: Combining Equation (13) with Equations ( 9)-( 12), it can be obtained that: in which it is quite easy to find that: Currently, many researchers [28,29] have concluded that the hardening laws for effective stress σ and effective plastic strain ε p were obtained from experimental data.The following power hardening laws [34] are adopted: In this study, the Langmuir function with three unknown parameters is introduced to characterize the hardening behavior of materials, and it can be defined as: where α, β, and n are the material parameters that can be obtained through nonlinear fitting of the effective stress-plastic strain curves.

Damage Initiation Criteria
In order to predict the initiation of various damage modes within a single layer and to evaluate the effective stress state during the loading process, this work uses the maximum stress criterion and the matrix failure criterion in the LaRC05 criterion to determine the damage initiation of fiber and matrix tension, respectively.
(1) Failure of the fiber tension (F ft ≥ 1, σ 11 > 0): where X T represents the tensile strength in the fiber direction of unidirectional FRP laminates.
(2) Failure of the matrix tension (F mt ≥ 1, σ n (θ T ) > 0): Matrix failure, also known as inter-fiber failure, refers to the occurrence of initial fracture surfaces in the matrix of unidirectional FRP laminates under transverse stress and in-plane shear stress.The subsequent damage and failure behavior of the matrix will propagate along the fracture surfaces parallel to the fibers.As shown in Figure 5, the various stress components acting on the fracture surface of the matrix are given, where 1-2-3 represents the material coordinate system of the unidirectional lamina, L-N-T represents the local coordinate system of the fracture surface, and the 1-axis coincides with the L-axis.θT is the rotation angle between the normal direction of the potential fracture surface of the matrix and the 2-axis of the material coordinate system, and the range of rotation angle θT is [−90 • , 90 • ]. σNN(θT) represents the normal stress on the fracture surface of the matrix, while σNT(θT) and σLN(θT) represent the shear stresses perpendicular and parallel to the fiber direction on the fracture surface of the matrix, respectively.Pinho et al. [35,36] established a criterion to determine the initial damage of the matrix in unidirectional FRPs laminates, based on the Mohr-Coulomb fracture theory, using stress components in the local coordinate system of the initial fracture surface of the matrix.
In Equation (19), ST and SL represent the shear strengths perpendicular to the fiber direction and along the fiber direction, respectively, on the fracture surfaces of the matrix; ηT and ηL are the transverse shear and longitudinal shear friction coefficient of the matrix, respectively; YT is the transverse tensile strength of the unidirectional FRPs laminate; and the relationship among ST, SL, ηT, ηL, and YC can be given by: In Equation (20), YC represents the transverse compressive strength of the unidirectional FRPs laminate; θC denotes the angle between the fracture surface of the matrix and the 3-axis in the material coordinate system, as failure occurs under transverse compression.For CFRP composite laminates, researchers such as Puck et al. [37] have experimentally determined its range of values to be [51°, 55°].
The following relationship exists between the stress on the fracture surface and the Pinho et al. [35,36] established a criterion to determine the initial damage of the matrix in unidirectional FRPs laminates, based on the Mohr-Coulomb fracture theory, using stress components in the local coordinate system of the initial fracture surface of the matrix.
In Equation ( 19), S T and S L represent the shear strengths perpendicular to the fiber direction and along the fiber direction, respectively, on the fracture surfaces of the matrix; η T and η L are the transverse shear and longitudinal shear friction coefficient of the matrix, respectively; Y T is the transverse tensile strength of the unidirectional FRPs laminate; and the relationship among S T , S L , η T , η L , and Y C can be given by: In Equation (20), Y C represents the transverse compressive strength of the unidirectional FRPs laminate; θ C denotes the angle between the fracture surface of the matrix and the 3-axis in the material coordinate system, as failure occurs under transverse compression.For CFRP composite laminates, researchers such as Puck et al. [37] have experimentally determined its range of values to be [51 The following relationship exists between the stress on the fracture surface and the initial fracture surface angle θ T of the matrix: (21) where:

Damage Evolution Law
When the internal stress of the material satisfies the failure criteria, the further increase in the effective stress will lead to the evolution of damage.During this process, strain energy is continuously released and the material exhibits localized softening phenomena, including the degradation of mechanical properties and a decrease in load-bearing capacity.By observing the experimental stress-strain curve of unidirectional T700/PEEK laminates, the curve segment concerning the damage evolution shows a sudden drop and degradation after reaching the ultimate tensile strength.In this study, a linear progressive degradation model is adopted to describe the damage evolution of the material [38,39].
where d t I represents the damage parameter at the current incremental step; ε is the equivalent strain in the composite ply and the strain values ε 0 I and ε f inal I are the equivalent strains corresponding to the initiation failure and final failure, respectively.When the material is undamaged, d t I = 0; when the material is completely failed, d t I = 1.In order to prevent the singularity of the material's stiffness matrix due to stiffness reduction to zero during the finite element calculation process, d * I = 0.999 is set.The material strain is the direct cause of damage evolution, so when the strain energy release rate of the material is equal to its fracture toughness, the material has experienced complete failure.As the material model behaves softer, the damage of the material exhibits localizing characteristics as the dissipated energy decreases upon mesh refinement [40].In order to avoid the effect of mesh size on the simulation results, Bazant's crack band theory [41] is chosen in this work and l*, which represents the characteristic length of the element, is introduced for the numerical simulation analysis considering the linear progressive degradation.The energy released per unit volume g I,m can be determined by the specific energy release rate G I,c of the material.
The equivalent strains ε I of various damage modes at any given time, at the beginning of initial damage ε 0 I , and at complete failure ε f inal I , as well as the equivalent stresses σ 0 I at the beginning of initial damage and the corresponding material fracture toughness calculation formula, are as follows: (1) Fiber Tensile Failure: For the fiber damage mode, the shear stress on the fracture surface is generally small.In this model, the consideration of equivalent strain is directly replaced by the strain in the direction of the fiber.
(2) Matrix tensile failure: where ⟨•⟩ represents the Macaulay operator, where ⟨x⟩ = (x + |x|)/2 when x ∈ R. G IC , G IIC , and G IIIC refer to the fracture toughness corresponding to the opening-type I, sliding-type II, and tearing-type III cracks in the matrix, respectively.ε NN is the normal strain on the fracture surface of the matrix.G mt,c is the equivalent fracture toughness of the matrix under mixed mode.ε NT and ε LN denote the shear strains that are perpendicular and parallel to the fiber direction on the fracture surface of the matrix, respectively.The calculation formulas are as follows: where: Under the tensile loading, as F mt ≥ 1, the matrix damage occurs, which means that cracks occur in the matrix of the potential fracture surface.The mentioned damage can be simulated by reducing the traction component on the potential fracture surface, defined as follows: (32) where σ d NN , σ d LN , and σ d NT in Equation ( 32) represent the softened normal stress and shear stress on the potential fracture surface of the matrix, respectively.In this paper, we assume that the composite material is transversely isotropic and we defined σ TT and σ LT to be softened in the same manner as described by Equation (32).The softened stress in the global coordinate system can be given by: When the fiber reaches its ultimate tensile strength, its damage failure also affects the damage evolution of the matrix.Therefore, the stress reduction calculation formula caused by fiber tensile failure is defined as: where σ d ij and σ d ′ ij are the stresses after the softening of the material, corresponding to matrix fracture and fiber damage, respectively, in the global coordinate system.

Model Validation 4.1. Implementation of the Numerical Simulation
As shown in Figure 6, the virtual specimen was discretized using eight-node linear reduced integration (C3D8R) solid elements with a size of 2.5 mm × 2.5 mm × 0.125 mm and assumes that the interfaces between layers are completely bonded.The left end of the finite element model is constrained to restrict the displacement in the loading direction and the right end loading is controlled by a displacement applied to a reference point outside its face.Computational accuracy was set as double precision to reduce the accumulation error during simulation.The proposed 3D elastoplastic damage constitutive model was implemented in ABAQUS 2021, using a user-defined material subroutine (VUMAT).A detailed flow chart of the VUMAT is presented in Figure 7.To validate the model, offaxis tensile analysis tests were simulated by using ABAQUS 2021.As shown in Table 2, the material parameters for T700/PEEK unidirectional laminates are presented, among which the fiber tensile fracture energy, matrix fracture energy, in-plane shear strength, and transverse compressive strength are taken from the literature [30,42,43].The remaining material parameters are measured from the manufacturer's experiments and are provided by the manufacturer.
(1 ) σ are the stresses after the softening of the material, correspond matrix fracture and fiber damage, respectively, in the global coordinate system.

Implementation of the Numerical Simulation
As shown in Figure 6, the virtual specimen was discretized using eight-node reduced integration (C3D8R) solid elements with a size of 2.5 mm × 2.5 mm × 0.12 and assumes that the interfaces between layers are completely bonded.The left end finite element model is constrained to restrict the displacement in the loading dir and the right end loading is controlled by a displacement applied to a reference outside its face.Computational accuracy was set as double precision to reduce the mulation error during simulation.The proposed 3D elastoplastic damage consti model was implemented in ABAQUS 2021, using a user-defined material subr (VUMAT).A detailed flow chart of the VUMAT is presented in Figure 7.To valida model, off-axis tensile analysis tests were simulated by using ABAQUS 2021.As s in Table 2, the material parameters for T700/PEEK unidirectional laminates are pres among which the fiber tensile fracture energy, matrix fracture energy, in-plane strength, and transverse compressive strength are taken from the literature [30,42,43 remaining material parameters are measured from the manufacturer's experimen are provided by the manufacturer.

Determination of the Plasticity Parameters
The relationship between the effective stress and effective plastic strain for material is uniquely determined and is independent of loading conditions [44].Th the plastic parameters of the constitutive model in the plastic phase can be a through longitudinal tensile tests.As shown in Figure 8, the composite laminate jected to a longitudinal tensile load, then a state of plane stress can be assumed.
The stress components are transformed along the principal material axes and pressed in terms of axial stress (σx), as follows:

Determination of the Plasticity Parameters
The relationship between the effective stress and effective plastic strain for a given material is uniquely determined and is independent of loading conditions [44].Therefore, the plastic parameters of the constitutive model in the plastic phase can be attained through longitudinal tensile tests.As shown in Figure 8, the composite laminate is subjected to a longitudinal tensile load, then a state of plane stress can be assumed.
The stress components are transformed along the principal material axes and are expressed in terms of axial stress (σ x ), as follows: where θ represents the angle between the tensile loading direction and the fiber direction, and the effective stress σ can be obtained by combining Equations ( 35)- (37) with Equation ( 8): σ = h(θ)σ x (38) and the transformation function h(θ) can be given by: Assuming that the studied material exhibits significant plastic deformation an deformation.In terms of the experimentally measured stress-strain data, the incr plastic strain p x dε can be obtained by subtracting the incremental elastic strain dε the incremental strain dεx, and it can be expressed as: dε dε dε The incremental elastic strain e x dε can be computed from: The elastic modulus Ex at different off-axis angles can be experimentally achi the stress-strain curve, and it also can be theoretically determined by the elastic m E1 and E2, off-axis angel θ, as well as the shear modulus G12 and Poisson 's ratio theoretical calculation formula can be given by: As shown in Figure 9, a comparison concerning the elastic modulus in loadin tion is made between the experimental data and calculation value obtained by E (42).It can be observed that there is a comparatively good agreement between th imental and theoretical values, which confirm the validity of the experiment in th Assuming that the studied material exhibits significant plastic deformation and small deformation.In terms of the experimentally measured stress-strain data, the incremental plastic strain dε p x can be obtained by subtracting the incremental elastic strain dε e x from the incremental strain dε x , and it can be expressed as: The incremental elastic strain dε e x can be computed from: The elastic modulus E x at different off-axis angles can be experimentally achieved by the stress-strain curve, and it also can be theoretically determined by the elastic modulus E 1 and E 2 , off-axis angel θ, as well as the shear modulus G 12 and Poisson 's ratio ν 12 .The theoretical calculation formula can be given by: As shown in Figure 9, a comparison concerning the elastic modulus in loading direction is made between the experimental data and calculation value obtained by Equation (42).It can be observed that there is a comparatively good agreement between the experimental and theoretical values, which confirm the validity of the experiment in this work.According to the coordinate transformation law of plastic strain components crement of plastic strain in the uniaxial loading direction can be expressed as [45]  Based on Equations ( 8), ( 10), ( 15) and ( 35)-( 37), the effective plastic strain can plified as: p ( ) For the determination of plasticity parameters a and b in Equation ( 8), many r ers adopt artificial methods, like trial and error, to obtain them.In this study, an o function regarding the plasticity parameters is provided to make each curve of e stress-plastic strain associated with different off-axis angles coincide with one curve.The plastic parameters a and b can be obtained by minimizing the objective f in Equation ( 45) by adopting the quasi-Newton algorithm.This aims to further e the consistency of the effective stress-plastic strain curves at different off-axis ang expressions are given as follows: 15 , 30 , 45 , 60 , 75 , 90 , where p ( ) θ σ ε is the effective plastic strain at different off-axis angles.
The plastic parameters were then calculated by using the quasi-Newton al programmed in MATLAB R2021b.The initial and final values of the effective plast were set to 10 −6 and 0.01, respectively, with a step size of 10 −6 .The initial value plastic parameters were set to 0.5, and the optimization results of the objective f are provided in Table 3.It can be observed that the optimal objective function v the two-parameter function model is 18,491.6292,whereas the optimal objective f value for the single-parameter function model is 31,926.5522.This indicates that posed approach is more suitable for characterizing the plastic behavior of CFR thermoplastic resins.According to the coordinate transformation law of plastic strain components, the increment of plastic strain in the uniaxial loading direction can be expressed as [45]: Based on Equations ( 8), ( 10), ( 15) and ( 35)- (37), the effective plastic strain can be simplified as: For the determination of plasticity parameters a and b in Equation ( 8), many researchers adopt artificial methods, like trial and error, to obtain them.In this study, an objective function regarding the plasticity parameters is provided to make each curve of effective stress-plastic strain associated with different off-axis angles coincide with one master curve.The plastic parameters a and b can be obtained by minimizing the objective function in Equation ( 45) by adopting the quasi-Newton algorithm.This aims to further enhance the consistency of the effective stress-plastic strain curves at different off-axis angles.The expressions are given as follows: where σ θ (ε p ) is the effective plastic strain at different off-axis angles.
The plastic parameters were then calculated by using the quasi-Newton algorithm programmed in MATLAB R2021b.The initial and final values of the effective plastic strain were set to 10 −6 and 0.01, respectively, with a step size of 10 −6 .The initial values of the plastic parameters were set to 0.5, and the optimization results of the objective function are provided in Table 3.It can be observed that the optimal objective function value for the two-parameter function model is 18,491.6292,whereas the optimal objective function value for the single-parameter function model is 31,926.5522.This indicates that the proposed approach is more suitable for characterizing the plastic behavior of CFRPs with thermoplastic resins.The power function of plastic hardening was regarded as the existing and generally adopted approach to characterizing the hardening law.In this work, the effective stress-plastic strain relationships are fitted by using the Langmuir function.According to Equations ( 38) and ( 44), the axial stress-plastic strain curves at different off-axis angles were transformed into effective stress-plastic strain curves, which were analyzed and fitted with a main curve by two-parameter Langmuir function (2 para_Lan) model, as shown in Figure 10.The processed data points show good consistency, and the parameters of the fitted hardening function are as follows: α = 200.01,β = 87.61,and n = 0.28.The power function of plastic hardening was regarded as the existing and generally adopted approach to characterizing the hardening law.In this work, the effective stressplastic strain relationships are fitted by using the Langmuir function.According to Equations (38) and ( 44), the axial stress-plastic strain curves at different off-axis angles were transformed into effective stress-plastic strain curves, which were analyzed and fitted with a main curve by two-parameter Langmuir function (2 para_Lan) model, as shown in Figure 10.The processed data points show good consistency, and the parameters of the fitted hardening function are as follows: α = 200.01,β = 87.61,and n = 0.28.

Prediction of the Results and Verification
As shown in Figure 11, the failure envelope curve, predicted by the present failure criterion, is in good agreement with the experiment data.Under the condition of the same loaded area of the specimen, the tensile strength decreases rapidly within the range of 0° to 15° of off-axis angle.However, as the off-axis angle further increases, the decreasing trend of the tensile strength gradually becomes gentler.The cloud map in Figure 12 shows the amount of matrix damage variables at the ultimate load and end of the analysis of the specimen, which was quantified by the variable of matrix damage dmt.As the dmt is 0, it indicates that the matrix has not been damaged, while when the dmt is 0.999, it indicates that the matrix has completely failed and has lost its bearing capacity.As is depicted in this cloud map, the matrix first suffers damage in the gage section of the specimen.With the further increase in the displacement load, the damage of the matrix begins to evolve towards the undamaged areas, and the evolution direction is consistent with the fiber direction.At the end of the analysis, the matrix damage has completely penetrated the crosssection of the specimen, directly leading to the complete failure and loss of load-bearing capacity of the specimen.

Prediction of the Results and Verification
As shown in Figure 11, the failure envelope curve, predicted by the present failure criterion, is in good agreement with the experiment data.Under the condition of the same loaded area of the specimen, the tensile strength decreases rapidly within the range of 0 • to 15 • of off-axis angle.However, as the off-axis angle further increases, the decreasing trend of the tensile strength gradually becomes gentler.The cloud map in Figure 12 shows the amount of matrix damage variables at the ultimate load and end of the analysis of the specimen, which was quantified by the variable of matrix damage d mt .As the d mt is 0, it indicates that the matrix has not been damaged, while when the d mt is 0.999, it indicates that the matrix has completely failed and has lost its bearing capacity.As is depicted in this cloud map, the matrix first suffers damage in the gage section of the specimen.With the further increase in the displacement load, the damage of the matrix begins to evolve towards the undamaged areas, and the evolution direction is consistent with the fiber direction.At the end of the analysis, the matrix damage has completely penetrated the cross-section of the specimen, directly leading to the complete failure and loss of load-bearing capacity of the specimen.The stress-strain curves obtained from numerical simulations and experiments for the 15° and 90° layups are compared in Figure 13.It can be observed that the proposed model in this study can effectively predict the mechanical behavior of the unidirectional T700/PEEK composite materials.Figure 14 shows the plastic strain components-total axial strain relations under different off-axis angles.Due to the assumption that there is no plastic deformation in the fiber direction, 1 1 p ε is always zero.The plastic strain components The stress-strain curves obtained from numerical simulations and experiments for the 15° and 90° layups are compared in Figure 13.It can be observed that the proposed model in this study can effectively predict the mechanical behavior of the unidirectional T700/PEEK composite materials.Figure 14 shows the plastic strain components-total axial strain relations under different off-axis angles.Due to the assumption that there is no plastic deformation in the fiber direction, 1 1 p ε is always zero.The plastic strain components The stress-strain curves obtained from numerical simulations and experiments for the 15 • and 90 • layups are compared in Figure 13.It can be observed that the proposed model in this study can effectively predict the mechanical behavior of the unidirectional T700/PEEK composite materials.Figure 14 shows the plastic strain components-total axial strain relations under different off-axis angles.Due to the assumption that there is no plastic deformation in the fiber direction, ε p 11 is always zero.The plastic strain components in the 1-3 and 2-3 plane, ε p 13 and ε p components which can be almost negligible.The curves of the plastic strain components, ε p 13 and ε p 23 , versus the total axial strain ε both like a straight line almost perpendicular to the vertical axis, and the two lines almost coincide with each other.It can be observed that the plastic strain components ε p 12 plays a major contribution on the plastic deformation for the off-axis angles 15 • .In particular, it can still be observed that the in-plane shear plastic strains generated are approximately 9 times and 12 times larger in the 2 and 3 directions, respectively.This indicates that plastic deformation is predominantly governed by the in-plane shear strain.However, as the off-axis angle increases further, accompanied by the weakening of in-plane plastic shear effects, the material's nonlinear behavior also gradually weakens.When the off-axis angle increases to 90 • , the in-plane plastic shear strain decreases to zero, and the plastic behavior of the material is mainly contributed to by the plastic strains in the 2 and 3 directions.plays a major contribution on the plastic deformation for the off-axis angles 15°.In particular, it can still be observed that the inplane shear plastic strains generated are approximately 9 times and 12 times larger in the 2 and 3 directions, respectively.This indicates that plastic deformation is predominantly governed by the in-plane shear strain.However, as the off-axis angle increases further, accompanied by the weakening of in-plane plastic shear effects, the material's nonlinear behavior also gradually weakens.When the off-axis angle increases to 90°, the in-plane plastic shear strain decreases to zero, and the plastic behavior of the material is mainly contributed to by the plastic strains in the 2 and 3 directions.plays a major contribution on the plastic deformation for the off-axis angles 15°.In particular, it can still be observed that the inplane shear plastic strains generated are approximately 9 times and 12 times larger in the 2 and 3 directions, respectively.This indicates that plastic deformation is predominantly governed by the in-plane shear strain.However, as the off-axis angle increases further, accompanied by the weakening of in-plane plastic shear effects, the material's nonlinear behavior also gradually weakens.When the off-axis angle increases to 90°, the in-plane plastic shear strain decreases to zero, and the plastic behavior of the material is mainly contributed to by the plastic strains in the 2 and 3 directions.(2) A two-parameter 3D plastic potential was developed by incorporating both deviatoric and dilatation deformations.The plastic parameters in the proposed plastic potential are determined using the quasi-Newton algorithm to optimize an objective function related to these plastic parameters.The Langmuir function was initially employed to fit the effective stress-strain relationships of the composites for hardening behavior, providing a precise prediction of the nonlinear mechanical response in CFRPs with thermoplastic resins.(3) The two-parameter 3D plastic potential was integrated with the Langmuir function of plastic hardening to model the constitutive behavior for unidirectional T700/PEEK Laminates.Based on the proposed constitutive model, a user-defined subroutine (VUMAT), based on ABAQUS, has been developed.This subroutine incorporates the maximum stress criterion for fiber failure and the LaRC05 criterion for matrix failure, allowing for the simulation of 3D elastoplastic damage behavior during tensile loading tests.This work constructs a two-parameter 3D elastoplastic damage constitutive model considering the material plasticity, damage evolution, and dilatational deformation.(4) Applying the 3D elastoplastic damage constitutive model to simulate the off-axis tensile behavior of unidirectional T700/PEEK laminates, the predicted stress-strain curves align closely with experimental data.This suggests that the model effectively captures the plastic behavior of unidirectional T700/PEEK thermoplastic laminates.The analysis of plastic strain in various directions indicates that shear plastic strain in the 1-2 plane governs plastic deformation.Nevertheless, this effect diminishes as the off-axis angle increases.The comparison between the numerically simulated and experimentally measured tensile strengths under identical loading planes reveals a gradual decrease in ultimate tensile strength with an increasing off-axis angle.

Conclusions
This work has investigated the nonlinear behavior of unidirectional thermoplastic composites under off-axis tensile loading.An anisotropic constitutive model for the nonlinear kinematic hardening elastoplastic damage of woven thermoplastic composites is still needed, including the asymmetry of tension and compression under cyclic loading conditions.

Figure 2 .
Figure 2. Experimental equipment for the off-axis tensile testing: (a) Overall diagram of the testing machine.(b) Specimen and extensometer.

Figure 2 .
Figure 2. Experimental equipment for the off-axis tensile testing: (a) Overall diagram of the testing machine.(b) Specimen and extensometer.
represents the plastic gradient, which describes the direction of the plastic strain increment component dε p ij .Materials 2024, 17, 3317 9 of 22

Materials 2024 , 23 Figure 5 .
Figure 5.The stress state of a single layer in the composite materials and the stress components on potential fracture surfaces.

Figure 5 .
Figure 5.The stress state of a single layer in the composite materials and the stress components on potential fracture surfaces.

Materials 2024 ,Figure 9 .
Figure 9.Comparison of the elastic modulus in loading direction between the theoretical an imental values under different off-axis angles.

Figure 9 .
Figure 9.Comparison of the elastic modulus in loading direction between the theoretical and experimental values under different off-axis angles.

Figure 10 .
Figure 10.Fitting of the effective stress-plastic strain curve under off-axis tensile loading.

Figure 10 .
Figure 10.Fitting of the effective stress-plastic strain curve under off-axis tensile loading.

Figure 11 .Figure 12 .
Figure 11.Comparison chart of the experimental and predicted tensile strength values at different off-axis angles.

Figure 11 . 23 Figure 11 .Figure 12 .
Figure 11.Comparison chart of the experimental and predicted tensile strength values at different off-axis angles.
the minimum values among all the plastic strain components which can be almost negligible.The curves of the plastic strain compothe total axial strain ε both like a straight line almost perpendicular to the vertical axis, and the two lines almost coincide with each other.It can be observed that the plastic strain components 1 2 p ε

Figure 13 .Figure 14 .
Figure 13.Comparison of the axial stress-strain curves between the experimental test and different prediction approaches under two off-axis angles: (a) 15°; (b) 90°.

( 1 )
By conducting off-axis tensile tests on seven types of samples, namely [0]8, [15]20, [30]20, [45]20, [60]20, [75]20, and [90]16, the corresponding axial stress-strain curves, elastic modulus, and tensile strengths were obtained.The fracture surface angles of the specimens for different off-axis angles were measured by using a high-precision digital angle gauge, and the difference between the measured fracture surface angle and the off-axis angle of each specimen was within −2°~3°.This difference is small and

Figure 13 .
Figure 13.Comparison of the axial stress-strain curves between the experimental test and different prediction approaches under two off-axis angles: (a) 15 • ; (b) 90 • .

Figure 13 .Figure 14 .
Figure 13.Comparison of the axial stress-strain curves between the experimental test and different prediction approaches under two off-axis angles: (a) 15°; (b) 90°.

( 1 )
By conducting off-axis tensile tests on seven types of samples, namely [0]8, [15]20, [30]20, [45]20, [60]20, [75]20, and [90]16, the corresponding axial stress-strain curves, elastic modulus, and tensile strengths were obtained.The fracture surface angles of the specimens for different off-axis angles were measured by using a high-precision digital angle gauge, and the difference between the measured fracture surface angle and the off-axis angle of each specimen was within −2°~3°.This difference is small and

Table 1 .
Fracture surface angles for the different layup methods.

Table 1 .
Fracture surface angles for the different layup methods.

Table 2 .
The material parameters of the T700/PEEK.

Table 2 .
The material parameters of the T700/PEEK.

Table 3 .
The calculation results of the optimization objective function between the single-parameter and two-parameter plastic potential.

Table 3 .
The calculation results of the optimization objective function between the single-parameter and two-parameter plastic potential.